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Abstract. We compute the effects of a compact flat universe on the angular 
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correlation function, the angular power spectrum, the circles-in-the-sky signature, and 
the covariance matrix of the spherical harmonics coefficients of the cosmic microwave 
background radiation using the full Boltzmann physics. Our analysis shows that the 
■ Wilkinson Microwave Anisotropy Probe (WMAP) three-year data are well compatible 

with the possibility that we live in a flat 3-torus with volume ~ 5 • lO'^Gpc'^. 



o 

. PACS numbers: 98.80.-k, 98.70.Vc, 98.80.Es 

^ ■ 1. Introduction 

O 

^ I At present, all data are consistent with, and in fact strongly support, the standard 
big bang model based on general relativity and the cosmological principle leading to 
^ I the general class of Friedmann-Lemaitre universes possessing the space-time structure 
I R X M^. Here R describes cosmic time and the three-dimensional comoving spatial 
^ ' section of constant curvature i^' = or i^' = ±1. Since the Einstein gravitational field 
equations are differential equations, they govern the local properties of space-time but 
^ I not the global geometry of the Universe at large, i. e. they do not determine the spatial 
I curvature and topology, and thus the shape of the Universe. 

Although a full quantum theory of gravity has not yet been estabhshed, it is 
commonly assumed that our Universe emerged from quantum fluctuations during 
the Planck era and that its spatial curvature and topology are invariants, i. e. have 
not changed ever since. The spatial curvature K can be inferred from astronomical 
observations of the total energy density etot of the Universe (at the present epoch) which 
is usually expressed in terms of the dimensionless density parameter Qtot = etot/ccrit- 
{K = for fitot = I, K = ±1 for fitot ^ 1- Qrit = 3Hy/{87iG), where 
Ho = 100/i(km/s)/Mpc denotes the Hubble constant.) 

Measurements of the temperature anisotropy of the cosmic microwave background 
(CMB) by NASA's satelhte Wilkinson Microwave Anisotropy Probe (WMAP) suggest 
(using the 3-year data combined with other observations) fltot — 1.0 [1]. Since this value 
for the density parameter is compatible with Euclidean geometry, we shall assume in 
this paper that our Universe is, indeed, fiat. (For a discussion of positively or negatively 
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curved universes, see e.g. [21 El S] and |S1 El E], respectively.) Having thus fixed the 
curvature, there remains the question of the topology of our Universe. 

It is a mathematical fact that fixing the curvature K does not determine uniquely 
the topology and thus the shape of the universe. Only if it is assumed that the spatial 
section is simply connected, it follows that the topology is (in the fiat case) given by the 
infinite Euclidean 3-space E^. Exactly this assumption is made in the popular ACDM 
model [1] which is composed of radiation (r), ordinary baryonic matter (b), cold dark 
matter (cdm) and has a cosmological constant A, i. e. Qtot = + + ^cdm + ^^a = 1- 

It is the purpose of this paper to demonstrate that the WMAP Syr data pj are 
well compatible with the possibility that we live in a "small Universe" having the 
shape of a fiat 3-torus whose fundamental domain is a cube with side length L ^ 4 
corresponding to a volume of ~ 5 ■ lO'^Gpc^. (L is given in units of the Hubble length 
Ljj = c/i^o — 4.26Gpc for h = 0.704.) For previous works on a toroidal universe, see 



P[91[I01[II1[I21[I31[T1[I51[I61[I71[T81[T91[2^ It turns out that 



the torus model describes the data much better than the best-fit ACDM model since it 
exhibits the suppression of the CMB anisotropy at large scales first observed by COBE 
[28]. 

Let us point out that Linde [29] has argued recently that a compact flat universe like 
the toroidal universe provides the simplest way to solve the problem of initial conditions 
for the low-scale inflation. While the quantum creation of a closed or an infinite open 
inflationary universe is exponentially suppressed, there is no such suppression for the 
toroidal universe. 

In [1] the WMAP team offers several sets of cosmological parameters depending 
on which cosmological data sets are taken into account. In this paper we use the 
cosmological parameters of the ACDM model of Table 2 in [Ij based on all astronomical 
observations, see their column "3 Year + ALL Mean", i.e. fib = 0.044, ficdm = 0.223, 

= 0.733, h = 0.704, Us = 0.947, r = 0.073. The distance to the surface of last 
scattering (SLS) is LgLs = Ar^Lj^ ~ 14.2Gpc where At] = tiq — rjsLs = 3.329 {rj is the 
conformal time). 

The temperature fluctuations 6T{n) of the CMB are calculated using the full 
Boltzmann physics of the coupled baryon-photon fluid, i. e. the ordinary and the 
integrated Sachs- Wolfe effect, the Doppler contribution, the Silk damping and the 
reionization are taken into consideration. The reionization is taken into account by 
modifying the ionization fraction Xe computed with recfast [30] 



by adding a term proportional to an erfc-function with parameters a and (3. (Here 
Yp = 0.24 denotes the Helium abundance.) We set a = 0.4 and (3 = 9 which leads to an 
optical depth of r = 0.073. 

The crucial difference between the toroidal and the simply connected ACDM 



(i/j., Uy, ^ due to the periodic boundary conditions imposed on the regular 




(1) 



model arises from the discrete spectrum of the vibrational modes {k 
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solutions of the Helmholtz wave equation for the 3-torus. (Note that a given wave 



number of representations of a e N by a sum of 3 squares of integers.) Since the 3-torus 
does not support wave lengths longer than L, there is an infrared cutoff k > 2tt/L (the 
mode = is not included being a pure gauge mode) which leads to the observed 
suppression of CMB power at large scales as we shall show now. In our calculations 
we use the exact mode spectrum for the first 10 000 eigenvalues, i. e. the first 5 503 602 
eigenmodes. The contribution of the higher modes is taken into account by a smooth 
remainder term as in the case of an infinite universe. 

2. Large scale temperature correlations and the angular power spectrum 

Let us consider the two-point temperature correlation function C{;d) = {6T{n) 6T{n')), 
h-n' = cosi), where (. . .) denotes an ensemble average over the primordial fluctuations, 
respectively, an ensemble average over the universal observers. In Fig. [1] we show as a 
function of the torus length L the integrated weighted temperature correlation difference 



between C^°'^'^^{'d) for the torus model and C°^^{i}) for four data sets derived from the 
WMAP 3yr data. The full and the dotted curves show I{L) derived from the ILC 3yr 
map and the TOH 3yr cleaned map [31] both with KpO mask, respectively, while the 
dashed and dotted-dashed curves refer to the ILC 3yr map and the TOH 3yr cleaned 
map both without mask. The horizontal lines show the corresponding results for the 
best-fit infinite ACDM model. It is seen from Fig. [T]that I{L) possesses a pronounced 
minimum at L = 3.86 for the maps with KpO mask, respectively at L = 4.35 for the 
maps without mask, which allows us to fix the best-fit torus model depending on the 
used data set. This shows that the WMAP-ILC 3yr map and the TOH 3yr cleaned 
map lead to the same best torus model. The uncertainty results mainly from the data 
near the galactic plane, i. e. by including them using no mask or by excluding them by 
applying the KpO mask. However, in both cases, the torus length L is close to L = 4, 
such that the main result is not affected by the foregrounds near the galactic plane. 
The values of I{L) obtained for the best-fit models are substantially lower than the 
corresponding values for the best-fit ACDM model demonstrating the suppressed CMB 
power of the toroidal universe. Obviously, if the toroidal universe gets too large (L > 9), 
the results are very similar to the infinite fiat universe. 

In Fig. [2] we show the mean temperature correlation function C (^9) for the torus 
model with the length L = 3.86 (full line) together with the la standard deviation (grey 
band). The la band is obtained from 10 000 sky map simulations of the considered torus 
model. For each realization the correlation function C(?9) is computed which in turn 
gives the standard deviation for this sample. This is compared to the WMAP-ILC 3yr 
(KpO) data (dotted line) and the best-fit ACDM model (dashed line). As expected, the 
torus model is in much better agreement with the data than the ACDM model. Notice 
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k I contributes with a highly irregular multiplicity determined by the 
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Figure 1. The integrated weiglited temperature correlation difference I{L) is shown 
for the torus model in dependence on the length L and for the best-fit ACDM model 
(straight horizontal lines) as described in p] incorporating the cosmological parameters 
of the "WMAPSyr + ALL" model. The models are compared with four data sets 
derived from the WMAP Syr data. 
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Figure 2. The temperature correlation function C(i?) is shown for the torus model 
with the length L ~ 3.86 and for the best-fit ACDM model. The grey band shows the 
1(7 deviation computed from 10000 simulations of the L — 3.86 torus model. 
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Figure 3. The angular power spectrum ST^ is shown for the torus model with the 
length L = 3.86 with the cosmic variance in comparison with the WMAP 3yr data. 

that the difference C°^^{i)) — (7™°dci^^^ between the data and the torus model is almost 
always a factor two smaller than the difference between the data and the ACDM model. 
In Fig. [3] we compare the power spectrum 6T^ = l{l + \)Ci/{2-k) (including cosmic 
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variance) for the same torus model with the WMAP Syr power spectrum [T], which is 
produced by combining the maximum hkehhood estimated spectrum from / = 2 — 10 
with the pseudo-C; based cross-power spectra for / > 10. This WMAP Syr power 
spectrum matches the power obtained from the previous sky maps without using any 
masks. Since we are mainly concerned with the large scale structure, we present only 
the low multipoles with 2 < / < 20. We have checked, however, that the power spectrum 
reproduces correctly the higher acoustic peak structure. The overall normalization is 
the only free parameter and has been determined by adjusting the power spectrum in 
the range 20 < / < 100 to the WMAP data. 

In Fig. m and Fig. [5] a similar comparison is carried out for the torus model with 
L = 4.S5 which now refers to the WMAP-ILC Syr data without mask. Again a 
remarkable agreement with the data is found. 
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Figure 4. The temperature correlation function C{-d) is shown for the torus model 
with the length L ~ 4.35 and for the best-fit ACDM model. The grey band shows the 
Icr deviation computed from 10000 simulations of the L — 4.35 torus model. 
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Figure 5. The angular power spectrum ST^ is shown for the torus model with the 
length L — 4.35 with the cosmic variance in comparison with the WMAP 3yr data. 
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3. The circles-in-the-sky signature 

Up to now we have compared the statistical properties for the CMB of a torus model 
with the observations. There are, however, two quantities, which can reveal the topology 
of the Universe more directly. Besides a signature in the covariance matrix, to which we 
devote the next section, this is the circles-in-the-sky signature [32] . On each circle pair on 
the CMB sky which is connected by the periodicity condition of the topological cell, there 
should be similar temperature fluctuations due to the Sachs- Wolfe effect. Expanding the 
temperature 5Tj(0) along a circle in a Fourier series 5Tj(0) = -^«ni^™"^' ^ ^ 27r, 
allows one to define the m-weighted circle signature for two circles i and j having a 
radius a as [22j 

2 V mT- T* p^™^ 

E.„"(|r„„P + |r,,„p) ■ 

The angle /3 takes a possible shift between the two circles into account. In [22] the 
WMAP lyr data are analysed and no topological signature is found. In order to 
investigate the question whether this negative result already rules out a torus cell as 
advertised in the previous section, we produce sky-map simulations for a torus with 
length L = 4 using the full Boltzmann physics and taking all eigenmodes up to the 
wavenumber k = 384.75 into account, i.e. in total 61,556,892 different wavevectors k. 
The contributions of the modes are computed up to /max = 1000. A torus model with 
L = 4 possesses 3 circle pairs with radius a ~ 53.08° and 6 circle pairs with radius 
a ~ 31.84°. The sky maps are computed in a HEALPix resolution iVgide = 512 and 
are then smoothed to a resolution of 0.5° and finally downgraded to Nsidc = 128. The 
maximum S{a) := maxijp Sij{a, P) of the circle signature ([3]) reveals the torus cell 
clearly as shown in Fig.[6l where the arrows indicate the peaks corresponding to the two 
radii. Now we disturb the sky maps by modifying the spherical expansion coefficients 

(4) 

V 1 + r ^ ^ 

where / is a constant factor and p{m) is a permutation of m G [— /, /]. With increasing / 
an increasingly disturbed sky map is generated without altering the statistical properties 
since only the I subspace is permuted. We found that the topological peaks in the circle 
signature S{a) vanish if the sky maps are modified by 50;uK on the average (see Fig.E]). 
Thus the crucial question emerges: what is the accuracy of the present CMB sky maps, 
i. e. is the CMB accuracy better than 50/iK? In the latter case the analysis carried out 
in [22] would rule out the torus model. The ILC map provided by the WMAP team does 
not specify the accuracy, and its usage for cosmological analyses is not recommended. 

Let us at first discuss the detector noise which is given by the WMAP 
team [33] for the sky maps for the 8 differencing assemblies x with x G 
{Qi, Q2, Vi, V2, Wi, W2, VF3, W4}. We construct from the 8 foreground reduced sky maps 
5Tx a mean value per pixel ST{i) := | J2x ^"^xif) where i denotes the pixel number with 
i G {0, 1, .., 12A^g^jjg — 1}. In order to carry out the average, the differencing assemblies 
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Figure 6. The circles-in-the-sky signature S{a) is sliown in dependence on the radius 
a of the back-to-back circles for a sky-map simulation of a torus with L — A (red) 
showing peaks at the corresponding circle radii a ~ 31.84° and a ~ 53.08° (indicated 
by arrows). Disturbing this sky map by 50/iK destroys this signature (blue). The 
result for the ILC-3yr map of the WMAP team is displayed as a green curve which 
possesses generally larger values than the simulation. 
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Table 1 . The average scatter ([6]) of the pixels outside the Kp2 mask is given in /iK for 
the lyr and the 3yr data with respect to the 8 differencing assemblies. The standard 
mean deviation of o-(i) from (cr(«))^P^ is also given. 




X are smoothed to the common resolution of the Qi band. With this mean value one 
can define a root mean square deviation per pixel 



\^-7 



measuring the scatter around the mean value ST{i) of a pixel. From this root mean 
square deviation a{i) one can define an average over a given domain which we take as 
the area outside the Kp2 mask of the 3 year WMAP data excluding also point sources, 
i.e. we define 

:= , (6) 

where the sums run only over the pixels outside the Kp2 mask. The results are presented 
in Table [1] for three different A^side resolutions for the lyr and the 3yr sky maps. In the 
case of the highest resolution Aside = 512 a mean scatter of 74/iK and 43/iK around 
the averaged temperature value is obtained for the lyr and the 3yr maps, respectively. 



Do we live in a "small Universe"? 



8 



frequency 


fore 


A3yr 
fore 


3yr 


K 






34.1 


Ka 






25.9 


Q 


66.3 


28.9 


21.3 


V 


33.1 


6.2 


16.7 


w 


50.5 


5.3 


22.1 



Table 2. The mean amplitude of the foreground which is taken into account in the 
lyr and Syr data outside the Kp2 mask is Hsted in the 2nd and 3rd column. The 
mean difference between the lyr and 3yr data outside the Kp2 mask is given in the 
last column. The values are given in /iK. 

Downgrading the maps to lower values of A^side averages neighbouring pixels and thus 
reduces the detector noise. The lowest resolution for which a circle signature S{a) can 
be computed is N^iAe = 128. Even in this case a mean scatter of 47/xK and 27yuK is 
obtained. Since the negative result in [22] with respect to a topological detection is based 
on the one-year data with a mean scatter of (cr(i))^^^ = 47/iK outside the Kp2 mask, 
a possible existing signature can with this accuracy easily be masked as the simulation 
presented in Fig. [6] demonstrates. 

Until now the focus was on the domain outside the Kp2 mask in order to obtain 
an estimate for the accuracy of the temperature pixel values ignoring a problematic 
foreground contamination. Within the Kp2 mask which covers 15% of the sky, the 
foreground dominates the CMB signal and an additional uncertainty with respect to 
the foreground removal arises. In the cleaned 3yr maps "substantial errors (> 30yuK) 
remain inside the Kp2 cut due to the limitations in the template model" ^3]. This 
residual foreground contamination is also the likely reason why we obtain different 
favoured lengths L of the torus cell depending on using the KpO mask or not as discussed 
in the previous section. An estimate for the uncertainty in the foreground removal 
procedure is obtained by taking the difference between the primary maps 5T(i)pri and 
the foreground reduced maps 5T(i)forered per frequency band provided by the WMAP 
team, i.e. by defining the average difference A^J^.^ := {5T{iYj^l — {i)]^^^^^^^^'^ and 
A^j;^ := {5T{ifj;^ - 5T(i)^j;^^^d)Kp2 outside the Kp2 mask. One obtains from Table [2] 
large differences in the component which is considered as foreground from the lyr to 
the 3yr data. The signal which is considered as foreground is much smaller in the 3yr 
data than in the lyr data. This change can be compared with the average difference 
A3y^ := (5T(i)p^[ — (5T(i)p^[)^P^ between the lyr and 3yr data which is also given in 
Table O An order of magnitude estimate of the foreground uncertainty may thus lie in 
the range 20 . . . 30/iK. 

Because of these uncertainties we conclude that the analysis in [22j does not rule 
out a toroidal universe. 
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4. Probabilities from the covariance matrix 

In [25] the covariance matrix C'^' := {cbimaa'm') torus universe is compared with 

the WMAP data and a 95% confidence limit is found that the torus length is L > 4. 
(Here (. . .) denotes an ensemble average.) The covariance matrix for the infinite model 
^im ~ Ci5ui5mm' being diagonal contrasts to the corresponding one of an anisotropic 
multiply connected space having non-vanishing non-diagonal elements. These can thus 
serve as a fingerprint for an anisotropic topology of the universe. However, this statement 
refers to the ensemble average. For an individual model non-diagonal elements are 
observed in both cases. In Fig. [7| the average S := ( (Cj'jT') ) over 

models of the non-diagonal elements is shown for several torus models and the infinite 
ACDM model. The grey band shows the average deviation from the mean value for the 
simulations belonging to the L = 3.86 torus. Here and in the following, the covariance 
matrix is truncated at Zmax = 30. Numerical tests based on sky maps of the torus 
model show that for side lengths around L ~ 4 the modes above / ~ 20 do not improve 
the likelihood to detect a torus topology. Therefore, the restriction to / < 30 in the 
covariance matrix is justified. For the infinite model, the average decays as 
leading asymptotically to vanishing non-diagonal elements. Fig. [7] shows that the torus 
models display indeed a deviation from the 1 / \/N decay leading to non- vanishing non- 
diagonal elements in the covariance matrix. To reveal this behaviour, an average over a 
large number of models is necessary, but we have observations only from one Universe. 
In addition, this is aggravated by the large average deviation as shown in Fig. [71 A 
closer inspection reveals that the non-vanishing non-diagonal elements in the case of the 
torus models are due to their real parts, whereas the imaginary parts vanish according 
to the 1 / Vn decay. 




Figure 7. The average S of the non-diagonal elements is shown in dependence on 
the ensemble size N for which the average is computed. The grey band displays the 
average deviation for the L = 3.86 torus ensemble. The 1/VN decay in the case of 
the infinite ACDM model is in this logarithmic plot revealed by the straight line with 
slope —1/2. The finite models deviate the stronger the smaller the torus length is in 
order to approach non- vanishing diagonal elements. 
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The covariance matrix can nevertheless be used to detect an anisotropic topology 
by computing the logarithm of the likelihood In £ = — ^(In det Css' + X^) + const., where 
s := /(/ + 1) + m and = Hiin J2s s' '^t^lJ/'^s'- The minimum is obtained by evaluating 
the sum ^, for all rotations of Og, i. e. for all possible orientations of the fundamental 
cell with respect to the CMB sky. The method of detecting a toroidal universe in this 
way is discussed in [26|l27j. In Fig. [8] we show In C for the covariance matrix of the torus 
model in dependence on the length L computed from three different sets of a/^ = cis- 
One set {as} is computed using a torus model with L = 3.86 (full curve) and provides 
the largest likelihood for the covariance matrix at L = 3.86 thus revealing the topology. 
No such maximum is observed using the obtained from the WMAP-ILC 3yr sky map 
(without mask) as shown by the dashed curve. Does this exclude a toroidal topology 
around L ~ 4? The answer depends on the accuracy of the obtained from the WMAP 
data. This accuracy in turn is determined by the accuracy of the ILC sky map, which 
we have discussed in the last section. To illustrate this point, we modify the phases of 
the as of the torus model with L = 3.86 by — > ■ 6""% where is a Gaussian random 
variable with unit variance. This blurs the phases by ~ 57°. The result is shown as 
the dotted curve in Fig. M and it is seen that this inaccuracy in the phases destroys the 
peak at L = 3.86. Such an inaccuracy in the phases would make the detection of an 
anisotropic topology impossible. The accuracy of the phases is discussed in [SU [35] and 
the WMAP team itself warns of using the ILC map for such a cosmological analysis. 
In [3l] it is shown that the Minimal Variance Method used in the construction of the 
ILC map leads to a wrong sign for ~ 40% of the coefficients a2o from 10^ Monte Carlo 
simulations, i.e. a phase shift of 180°. The assumption that the CMB is a Gaussian 
random field leads to phases which are uniformly distributed over the interval [0°, 360°]. 
In [35] it is shown that the phases for / < 10 are consistent with a uniform distribution 
but not the differences between the phases which should also be uniformly distributed. 
This should be taken as a warning that the phases might not be sufficiently accurately 
determined by the ILC map. 

Since the accuracy of the phases is unknown, we try to obtain an order of magnitude 
estimate by the following procedure. Assuming that the ILC map displays outside the 
kpO mask the genuine CMB signal, i. e. without any foregrounds, which is surely not 
true, and within the kpO mask a strongly disturbed signal which eradicates the CMB 
signal, we construct sky maps which have outside the kpO mask the corresponding part 
of a torus simulation and within the mask a map obtained from an infinite ACDM 
simulation. In this way we combine 950 torus sky maps with 100 infinite ACDM sky 
maps leading to 95 000 data sets of the spherical expansion coefficients aim- Comparing 
the phases of these sky maps with the phases of the corresponding torus simulations 
without mask reveals that the phase change A0 has a standard deviation o"(A0) ~ 63°. 
The large deviations which are of the same order as used in Fig. [8], have their roots in the 
strong oscillatory nature of the spherical harmonics which leads to many cancellations in 
the integration required for the computation of the a^m's. In Fig.[9K the distributions of 
the phase changes Acj) are shown for four different values of /, i. e. / = 5, 10, 20, and 30. 
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Figure 8. The logarithm of the hkehhood ln£ is shown for the covariance matrix Css' 
befonging to the torus model with length L. The is computed using the belonging 
to the torus model with L — 3.86 (full curve) and for the obtained from the WMAP- 
ILC Syr map (dashed curve). The sharp maximum in the full curve clearly "detects" 
the underlying torus model. However, modifying the phases destroys the maximum 
(dotted curve). 




The distributions are very similar and have a standard deviation cr(A0) independent of 
/ as revealed by Fig.[9}3. This demonstrates that the accuracy of the phases might be 
low enough to hide a signature of the topology. 

Since the true accuracy of the phases is unknown, one cannot draw a final conclusion 
whether the covariance matrix signature already rejects a toroidal topology for the 
Universe. Notice that the correlation function (^(-i?) and the angular power spectrum 
6T^ do not depend on these phases such that the very low power at angles {} > 60° is a 
much more robust observation pointing to a non-trivial topology. 




Figure 9. Panel a) shows the distribution of the error A(/) of the phases (in radian) 
by replacing the torus temperature fluctuations within the kpO mask by those of an 
infinite ACDM model for I = 5, 10, 20, and 30. The distribution with the highest 
peak at = belongs to ^ = 5 whereas the other distributions are indistinguishable. 
Panel b) shows that the standard deviation (j{^(j>) is cr(A0) ~ 1.1 independent of /. 
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The signature of the non-diagonal elements of the covariance matrix Cf^' is not only 
revealed in the spherical harmonics basis but also directly in the pixel space C(n, n'). 
This alternative is introduced in [36] and applied to compact hyperbolic spaces. In [37] 
a Bayesian analysis is carried out for the torus model using the pixel space correlation 
C{n,n') based on the COBE-DMR data. It is found that for e = L/(2Ar/) > 0.6, i.e. 
L > 4, the torus models are consistent with the standard ACDM model at the 2a level, 
and for some limited choices of orientations the fit to the COBE data is considerably 
better than that of the infinite model. Such a Bayesian analysis based on C{n, n') using 
the WMAP 3yr data is not carried out in this paper such that the advantage of the 
torus model over the ACDM model with respect to this measure remains open here. 

5. Summary 

Our analysis shows that the simplest non-trivial fiat topology, the cubic torus with a 
side length smaller than the diameter of the decoupling surface, is compatible with the 
WMAP 3-year data, and describes these data much better than the standard ACDM 
model. It is very intriguing that the toroidal topology provides also the simplest way to 
solve the problem of initial conditions for the low-scale inflation |29j . 

CMBFAST (www.cmbfast.org) and HEALPix (healpix.jpl.nasa.gov) [38] were used 
in this work. 
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